#include "hrv_timecompute.h"


HRV_TimeCompute::HRV_TimeCompute()
{

}

HRV_TimeCompute::TimeDomainHRV HRV_TimeCompute::timeDomainHRV(const QVector<double>& ibi, double win, double xx) {
    TimeDomainHRV output;

    // 计算基本统计量
    output.max  = *std::max_element(ibi.begin(),  ibi.end());
    output.min  = *std::min_element(ibi.begin(),  ibi.end());
    output.mean  = std::accumulate(ibi.begin(),  ibi.end(),  0.0) / ibi.size();
    output.median  = [&ibi]() {
        QVector<double> sortedIbi = ibi;
        std::sort(sortedIbi.begin(),  sortedIbi.end());
        int mid = sortedIbi.size()  / 2;
        return (sortedIbi.size()  % 2 == 0) ? (sortedIbi[mid - 1] + sortedIbi[mid]) / 2 : sortedIbi[mid];
    }();
    output.SDNN = [&ibi, &output]() {
        double variance = 0.0;
        for (double val : ibi) {
            variance += std::pow(val - output.mean,  2);
        }
        return std::sqrt(variance / ibi.size());
    }();
    output.SDANN = SDANN(ibi, win * 1000);
    auto [p, n] = pNNx(ibi, xx);
    output.NNx = n;
    output.pNNx = p;
    output.RMSSD = RMSSD(ibi);
    output.SDNNi = SDNNi(ibi, win * 1000);

    // 计算心率
    QVector<double> hr;
    for (double val : ibi) {
        hr.push_back(60.0  / (val / 1000.0));
    }
    output.meanHR  = std::accumulate(hr.begin(),  hr.end(),  0.0) / hr.size();
    output.sdHR  = [&hr, &output]() {
        double variance = 0.0;
        for (double val : hr) {
            variance += std::pow(val - output.meanHR,  2);
        }
        return std::sqrt(variance / hr.size());
    }();

    // 几何HRV计算
    int nbin = 32; // 默认分箱数
    output.HRVTi = hrvti(ibi, nbin);
    output.TINN = tinn(ibi, nbin);

    return output;
}

double HRV_TimeCompute::SDANN(const QVector<double>& ibi, double t) {
    QVector<double> means;
    int i1 = 0;
    double sum = 0.0;

    for (int i2 = 0; i2 < ibi.size();  ++i2) {
        sum += ibi[i2];
        if (sum >= t) {
            means.push_back(std::accumulate(ibi.begin()  + i1, ibi.begin()  + i2 + 1, 0.0) / (i2 - i1 + 1));
            i1 = i2 + 1;
            sum = 0.0;
        }
    }

    double mean = std::accumulate(means.begin(),  means.end(),  0.0) / means.size();
    double variance = 0.0;
    for (double val : means) {
        variance += std::pow(val - mean, 2);
    }
    return std::sqrt(variance / means.size());
}

double HRV_TimeCompute::SDNNi(const QVector<double>& ibi, double t) {
    QVector<double> stds;
    int i1 = 0;
    double sum = 0.0;

    for (int i2 = 0; i2 < ibi.size();  ++i2) {
        sum += ibi[i2];
        if (sum >= t) {
            double mean = std::accumulate(ibi.begin()  + i1, ibi.begin()  + i2 + 1, 0.0) / (i2 - i1 + 1);
            double variance = 0.0;
            for (int j = i1; j <= i2; ++j) {
                variance += std::pow(ibi[j] - mean, 2);
            }
            stds.push_back(std::sqrt(variance  / (i2 - i1 + 1)));
            i1 = i2 + 1;
            sum = 0.0;
        }
    }

    return std::accumulate(stds.begin(),  stds.end(),  0.0) / stds.size();
}

std::pair<double, double> HRV_TimeCompute::pNNx(const QVector<double>& ibi, double x) {
    QVector<double> differences;
    for (int i = 1; i < ibi.size();  ++i) {
        differences.push_back(std::abs(ibi[i]  - ibi[i - 1]));
    }

    int n = std::count_if(differences.begin(),  differences.end(),  [x](double val) { return val > x; });
    double p = (static_cast<double>(n) / differences.size())  * 100.0;
    return {p, static_cast<double>(n)};
}

double HRV_TimeCompute::RMSSD(const QVector<double>& ibi) {
    QVector<double> differences;
    for (int i = 1; i < ibi.size();  ++i) {
        differences.push_back(std::pow(ibi[i]  - ibi[i - 1], 2));
    }

    double sum = std::accumulate(differences.begin(),  differences.end(),  0.0);
    return std::sqrt(sum / differences.size());
}

double HRV_TimeCompute::hrvti(const QVector<double>& ibi, int nbin) {
    double minVal = *std::min_element(ibi.begin(),  ibi.end());
    double maxVal = *std::max_element(ibi.begin(),  ibi.end());
    double binWidth = (maxVal - minVal) / nbin;

    QVector<int> hist(nbin, 0);
    for (double val : ibi) {
        int bin = static_cast<int>((val - minVal) / binWidth);
        if (bin >= nbin) bin = nbin - 1;
        hist[bin]++;
    }

    int maxCount = *std::max_element(hist.begin(),  hist.end());
    return static_cast<double>(ibi.size())  / maxCount;
}

double HRV_TimeCompute::tinn(const QVector<double>& ibi, int nbin) {
    double minVal = *std::min_element(ibi.begin(),  ibi.end());
    double maxVal = *std::max_element(ibi.begin(),  ibi.end());
    double binWidth = (maxVal - minVal) / nbin;

    QVector<int> hist(nbin, 0);
    for (double val : ibi) {
        int bin = static_cast<int>((val - minVal) / binWidth);
        if (bin >= nbin) bin = nbin - 1;
        hist[bin]++;
    }

    int peakBin = std::distance(hist.begin(),  std::max_element(hist.begin(),  hist.end()));
    double minDiff = std::numeric_limits<double>::max();
    double tinnValue = 0.0;

    for (int m = peakBin - 1; m >= 0; --m) {
        for (int n = peakBin + 1; n < nbin; ++n) {
            double diff = 0.0;
            for (int i = m; i <= n; ++i) {
                double q = (i <= peakBin) ? (static_cast<double>(i - m) / (peakBin - m)) * hist[peakBin]
                                          : (static_cast<double>(n - i) / (n - peakBin)) * hist[peakBin];
                diff += std::pow(hist[i] - q, 2);
            }
            if (diff < minDiff) {
                minDiff = diff;
                tinnValue = std::abs((n - m) * binWidth);
            }
        }
    }

    return tinnValue;
}
